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The isothermal compressibility of a general crystal is analyzed within classical density functional 
theory. Our approach can be used for homogeneous and unstrained crystals containing an arbitrarily 
high density of local defects. We start by coarse-graining the microscopic particle density and 
then obtain the long wavelength limits of the correlation functions of elasticity theory and the 
thermodynamic derivatives. We explicitly show that the long wavelength limit of the microscopic 
density correlation function differs from the isothermal compressibility. It also cannot be obtained 
from the static structure factor measured in a scattering experiment. We apply our theory to crystals 
consisting of soft particles which can multiply occupy lattice sites (’cluster crystals’). The multiple 
occupancy results in a strong local disorder over an extended range of temperatures. We determine 
the cluster crystals’ isothermal compressibility, the fluctuations of the lattice occupation numbers 
and their correlation functions, and the dispersion relations. We also discuss their low-temperature 
phase diagram. 
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I. INTRODUCTION 

In crystals, where translational invariance is sponta¬ 
neously broken, strain enters as additional thermody¬ 
namic variable in the free energy, describing the distor¬ 
tion of the solid. The trace of the strain tensor is con¬ 
nected to the change in density. In particular, in ideal 
crystals, where all atoms can be unambiguously assigned 
to lattice sites and all lattice sites are occupied, density 
change is determined by the trace of the strain tensor. 
In real crystals, point defects like interstitials and vacan¬ 
cies are present, and density can change by both defor¬ 
mation of the solid (captured by the strain), and diffu¬ 
sion of defects 1 . Thus, the presence of defects opens the 
question how density and strain fluctuations are defined 
in real crystals. Here, no one-to-one mapping of atoms 
to lattice positions is possible. Therefore, the displace¬ 
ment field, whose symmetrized (in linear approximation) 
gradient gives the strain, cannot be obtained from the 
displacements of individual atoms from their lattice po¬ 
sitions. Only recently microscopic definitions of strain 
and density fluctuations in real crystals were derived 
from the statistical mechanical description of real crys¬ 
tals, overcoming this difficult}^. This work followed an 
earlier suggestion by Szamel and ErnstP 'A Preliminary 
Monte Carlo simulations and comparisons with older ap¬ 
proaches, including to amorphous solids, indicated the 
potential of the microscopic theorjEl. 

An intriguing finding of the microscopic approach of 
Ref. [2, concerns the coarse-grained density field Sn(r,t) 
which enters into the theory of crystal elasticity. Even 
for arbitrarily large wavelengths, particle density fluc¬ 
tuations with wavevectors close to all (finite) recipro¬ 
cal lattice vectors contribute to the coarse-grained den¬ 
sity field. In this contribution, we discuss this at first 
surprising finding within the framework of density func¬ 


tional theory. This theory allows us to properly link mi¬ 
croscopic and macroscopic density fluctuations in states 
with broken translational symmetry in order to parallel 
the coarse-graining of the free energy functional achieved 
previously for e.g. homogeneous liquid crystal^. 

Based on the microscopic definition of the coarse¬ 
grained variables of elasticity theory, we can ad¬ 
dress another intriguing question, originally raised by 

StillingeiEHn] 

Namely, whether the structure factor is 
an analytic function around zero wavevector and whether 
its small wavevector limit coincides with the compress¬ 
ibility? We find that due to the long-ranged displace¬ 
ment correlations, the small wavevector limit of the cor¬ 
relation function of the coarse-grained density field is 
non-analytic and depends on the direction relative to 
the crystal lattice. We derive these results from den¬ 
sity functional theory and can thus put them on a firm 
microscopic basis. Thus, we generalize earlier findings 
obtained within the harmonic crystal approximation 12 . 
Because of the non-analyticity, special care is required 
when discussing the thermodynamic limit. From studies 
on two-dimensional crystals it is known that defects en¬ 
ter the expression for the isothermal compressibility in 
a complicated fashion 1 ^. We generalize these results to 
crystals of arbitrary symmetry. Correcting the appendix 
of Ref. [2], we also derive relations between fluctuation 
functions and thermodynamic derivatives. These results 
suggest that the elastic constants of crystals with point 
defect:^ could be measured by microscopy techniques 
applied to colloidal crystal^. 

In order to t est th e theory, we apply it to so-called 
‘cluster crystals®!- 17 which consist of particles interact¬ 
ing with a soft-core repulsion. The softness of the po¬ 
tential allows for multiple occupancy of individual lat¬ 
tice sites by the particles and for fluctuations of the lat¬ 
tice sites occupation numbers. These fluctuations play 
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the role of mobile local defects. Indeed, the approxima¬ 
tion which considers these cluster crystals as ideal crys¬ 
tals (with a uniquely occupied lattice sites) is valid only 
at extremely low temperatures 1 ^, and the different crys¬ 
tal structures can only be understood by_ allowing for a 
distribution of site occupation number ::! 1 ! d 2 ° l. For these 
crystals, we will derive thermodynamic derivatives, cor¬ 
relation functions, and dispersion relations, which were 
not accessible previously, and we will discuss their low 
temperature phases. 

The paper is organized as follows: in Sect. II we first 
recall definitions and results from Ref. [2j and then, in 
Sect. Ill, we derive expressions for the fluctuations of 
displacement and density fields in an unstressed reference 
state. They are given by microscopic quantities defined in 
terms of the direct correlation function of the crystal. To 
facilitate application of these expressions, we also invert 
these relations considering two sets of independent fluc¬ 
tuations, coarse-grained density and displacement field or 
defect densit^ 21 ^ and displacement field. In Sect. IV we 
derive the thermodynamic free energy, including the ther¬ 
modynamic elastic susceptibilities, by coarse-graining the 
microscopic classical density functional. As the first step, 
we obtain the free energy functional containing the elas¬ 
tic fields, which reduces to the thermodynamic one for 
homogeneous fields. This is followed by the discussion 
of thermodynamic derivatives. In Sect. V we discuss the 
small wavevector limit of the coarse-grained density fluc¬ 
tuation function and show that it differs from the isother¬ 
mal compressibility n. We also discuss scattering func¬ 
tions and conclude that scattering experiments do not 
allow to measure the compressibility in a crystal, in con¬ 
trast to liquids and gase^2). Finally, in Sect. VI we apply 
our theory to cluster crystals. We show that a simple 
mean-field density functional leads to surprisingly accu¬ 
rate values of compressibilities and occupation number 
fluctuations. Details of some of the calculations are pre¬ 
sented in appendices. 


II. COARSE-GRAINED FIELDS 

Crystals exhibit spontaneously broken translational 
symmetry ( e.g. the average density is non-uniform) and 
this, via the Goldstone theorem, leads to long-ranged 
correlations. Specifically, the vector displacement field 
u(r, t ) possesses correlations which decay like the in¬ 
verse distance. In ideal crystals, one can use the familiar 
expression for the microscopic density of the displace¬ 
ment field Ui(t)5(r — Ri), involving the displacement 
Uj(f) = rj(f) — R, of the instantaneous position of the 
particle i , rj(f), from its lattice site R;. However, in real 
crystals, in which defects are present, this expression is 
invalid 3 . In order to find the microscopic definition for 
the displacement field u(r,f) and for the other fields of 
elasticity theory, an alternative approach was developed 
in Ref. [2]. 

Before we discuss the approach of Ref. [2] , we need to 


define precisely various fields used in the present paper. 
First, we have microscopic densities, i.e. quantities that 
are defined for and depend on an individual configuration 
of the iV-particle system. To distinguish these quantities 
we will always explicitly state that they depend on time t 
(like, e.g., in the standard definition of the displacement 
field mentioned in the previous paragraph.) Another ex¬ 
ample, which will be important in the following, is the 
microscopic particle density p(r, f); it will be precisely de¬ 
fined in Eq. ([T]) below. In crystals, in general the averages 
of microscopic quantities will change on the spatial scale 
of the crystalline cell. For example, the average density 
in a crystal, n(r) = (p(r,t)), is non-uniform, with large 
peaks near lattice sites’ positions. In contrast, the scalar 
density, denoted Sn(r) and the vector displacement field, 
<5u(r) used in the theory of elasticity vary only on much 
larger scales; here the S indicates a deviation from homo¬ 
geneous thermal equilibrium. Thus, one of the goals of 
Ref. [2] was to identify microscopic fields whose averages 
correspond to the fields of elasticity theory. In the rest 
of this paper we will call these fields microscopic coarse¬ 
grained fields. Also, in the rest of the paper we will re¬ 
fer to averages of microscopic quantities as macroscopic 
fields. Especially, second moments, viz. covariances and 
correlation functions, will be considered in the following 
and will be connected to thermodynamic derivatives. 


A. Microscopic particle density 

The concepts of generalized elasticity theorji 1 indi¬ 
cate that density fluctuations close to (all) reciprocal lat¬ 
tice vectors are long-rangecPH. Therefore, they all could 
contribute to coarse-grained fields. The microscopic ap¬ 
proach to find the displacement field in a real crystal 
starts from the particle density field p( r, t) which depends 
on the configuration of the IV-particle system (consid¬ 
ering, for simplicity, a one-component crystal of point 
particles interacting with a spherically symmetric pair- 
potential) 

N 

p(. r ,t) = j2 5 ( r w 

i= 1 

where r i(t) are the particle positions, N is the number 
of particles in the volume V ; later on we will use no to 
denote the average density, no = N/V. Spatial Fourier 
transformation gives fluctuations close to vectors g of the 
reciprocal lattice 

8p e (<l,t) =p(g + q,f) -n g IAS q ,o , (2) 

where 

N 

p(k,t) = / d\e-^ r p{v,t) = ^ e - ik r ‘W , 


( 3 ) 




and 


C. Relating the coarse-grained fields to 
microscopic density fluctuations 
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1 1 ^ 

n g = -(p( gj t)) = -^<e-^(*)). (4 ) 

i 

Here, the general wave vector k was divided up into recip¬ 
rocal lattice vector g and wave vector q, which lies within 
the first Brillouin zone; () brackets indicate canonical 
averaging at fixed temperature T (averages are time in¬ 
dependent due to time-translational invariance of equi¬ 
librium stately. n g are the Bragg-peak amplitudes 
(Debye-Waller factors) which serve as crystal order pa¬ 
rameters. They quantify the spontaneous breaking of the 
translational invariance (spatial homogeneity). 


B. Coarse-graining microscopic density 
fluctuations to elasticity fields 

In Ref. [2] the following representation was established 
for the microscopic density fluctuation in terms of micro¬ 
scopic coarse-grained density and displacement fields 

Sp g (q,t) = —in g g a 5u a (q, t) + , (5) 

n o 

with Greek indices denoting spatial directions; repeated 
indices are summed over (Einstein summation conven¬ 
tion is used). Equation ([5]) is the crucial relation linking 
the fields of macroscopic elasticity theory to the under¬ 
lying microscopic density fluctuations. It states that for 
wave vectors q within the first Brillouin zone, the four 
coarse-grained fields Sn(q, t) and du(q, t) determine the 
hydrodynamic contributions of the microscopic density 
field. This is valid even close to Bragg-peaks at arbi¬ 
trarily high reciprocal lattice vectors g. Equation ([5| 
was deduced considering the Zwanzig-Mori equations of 
motion of the microscopic density fluctuation^ 2 . In the 
present contribution, we support it by considerations of 
equilibrium correlations. 

In ideal crystals without defects the coarse-grained 
density and the divergence of the displacement field are 
proportional^. In real crystals, mass transport can arise 
from lattice distortions (described by the displacement 
field) but also from defect motion, which occurs diffu¬ 
sively over large distances. This additional hydrody¬ 
namic mode is called point defect density. It enters by 
the standard definition 1 : 

dc(q, t) = —<5n(q, t ) - in 0 q a Su a (q, t ) . (6) 

In Ref. [ 2 ] it is shown that Eqs. ([ 5 | and © predict the 
correct reversible dynamics of the defect density. Be¬ 
cause many situations require theoretical expressions at 
constant defect density®!, we will use Eq. © repeatedly 
in the following sections. 


Explicit expressions for the coarse-grained density and 
displacement fields can be derived by inverting Eq. ©• 
The inversion can be performed using the two summa¬ 
tions 


n 0 

Afo 



(7a) 

g 


] n g gp ■ 

(7b) 


and the relation )T) g |n g | 2 g = 0. The normalization 
constants are J\f 0 = J2 S \ n s\ 2 and AT a p = J2 S l n g| 2 £GS/3- 
Performing the sums over the reciprocal lattice in Eq. ([5]) 
leads to the microscopic coarse-grained density 

M q , t ) = ^Yl n s <W q ' *) > ( 8 ) 


and to the microscopic coarse-grained displacement field 
Su a (q, t ) = iM~p ^2 n *e 9P <5Pg( q , t) . (9) 

g 


These expressions could be evaluated using information 
obtained from computer simulations or from colloidal 
experimental 

Equations <[sj) and © express the coarse-grained fields 
in terms of microscopic particle density © . It is in¬ 
triguing that contributions from all finite lattice vectors 
g ^ 0 are present in the coarse-grained density. Even in 
the limit of vanishing wave vector, q —> 0, it is not suf¬ 
ficient to measure particle density fluctuations close to 
the center of the first Brillouin zone, in order to deter¬ 
mine the thermodynamic density field in crystals. Fluc¬ 
tuations from the regions around all lattice vectors con¬ 
tribute and describe how macroscopic strain fluctuations 
and defect density independently cause changes in the 
hydrodynamic particle density. 


III. RELATIONS INVOLVING CORRELATIONS 
OF THE COARSE-GRAINED FIELDS 

A. Correlation functions of the coarse-grained 
fields 

After recalling the relations between the fields of elas¬ 
ticity theory and microscopic fluctuation^, we turn now 
to the focus of our work, the correlations functions of 
the coarse-grained fields and the thermodynamic deriva¬ 
tives (including the isothermal compressibility). First, 
we will obtain the correlation functions from classical 
density functional theory (DFT)PH 2 ^. These correlation 
functions will then be analyzed in the homogeneous case 
to obtain the thermodynamic quantities. 
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Close to equilibrium, owing to the fluctuation dissipa¬ 
tion theorem, only equilibrium correlation functions are 
required in order to discuss the linear response to small 
external field^ 2 ®! In a homogeneous and unstrained crys¬ 
tal, the equilibrium correlation functions of the micro¬ 
scopic density fluctuations on the left hand side of Eq. ^ 
can be calculated within DFT. This enables us to obtain 
the correlation functions of the coarse-grained fields in 
Sects. ITTTAli . and IIII AH . The fundamental Ornstein- 
Zernike relation provides a connection between the den¬ 
sity correlations and the inverse density-density correla¬ 
tion matrix J gg / (q) 

Vk B T5ggn = y^(i5p g (q, t) 5pg> (q, t)) </ g ' g "(q) ■ (10) 

s' 

Here, the periodicity of the two-point density correlation 
functional was used which implies that only density fluc¬ 
tuations whose wavevectors differ by a vector of the re¬ 
ciprocal lattice are correlated. The (infinite-dimensional) 
Hermitian matrix J gg / is the double Fourier-transform of 
the second functional derivative of the free energy with 
respect to the macroscopic density, which includes as 
non-trivial part the direct correlation function c(ri,r 2 ). 


^gg'(q) — 

k B T 


( 11 ) 


d £ Vi/(I £ V2e ,;sri e^ iS '' r2 e 4qri2 


12) 

W r i) 


— c(ri,r 2 


The direct correlation f uncti on c(ri, r 2 ) is one of the cen¬ 
tral quantities of DFd^ 2 ^ and is obtained as second 
functional derivative of the excess free energy T ex with 
respect to the average density profile, fcgTc(r 1 ,r 2 ) = 

d ^ —^l)L Given an (approximate) expression for the 


<5n(ri)<5n(r2) 

free energy, J gg / can thus be taken as known. It con¬ 
stitutes the only input for the ensuing theory. As one 


consequence, in Sect. IV below only the quadratic ex¬ 
pression of the free energy functional will play a role and 
will be sufficient to evaluate the thermodynamic deriva¬ 
tives required for the elastic response. 


1. Including coarse-grained density 


It is now conceptually straightforward albeit somewhat 
tedious to derive the correlation functions of the coarse¬ 
grained fields in terms of expressions containing the di¬ 
rect correlation function. Using Eq. ([ 5 ]), one gets 


(5p* s {q.t)6pg'{q,t)) = 
n* s rig> (g a g' 0 {8u* a (q, t)5up( q, t)) 


( 12 ) 


,5n*(q, t)Sn( q, t ), 


ns 


+ i9a(du* a (q,t) 


(5?r(q, t) Sn*(q,t ) 


n 0 




n 0 


5u 0 {q,t))g'^j 


Inserting this into Eq. (10) and with the help of the two 
summations ([7]) and Eqs. ([8| and ([9]), one obtains the 


following set of equations 
(Sn*Sn) 


Vk R T = 


**( q)- 


0/3 = 


( Sn*Sn) 


Tin 


! Sn” 
v n 0 
8n 
K n 0 


-8u 0 )p f g(q), 

(13a) 

-Su s )X* 50 {q), 

(13b) 

* dn , 

l a ) v (q), 
n 0 

(13c) 

<7)K(q)' 

no 

(13d) 


Here, generalized (viz. ^-dependent) constants of elastic¬ 
ity, n(q),/z Q (q), and A a ^(q), appear. We will show that 
they enter into the equilibrium correlation functions of 
the coarse-grained fields and reduce to thermodynamic 
derivatives in the long-wavelength limit 2 . Using Eq. (Ill, 
the q-dependent constants of elasticity can be expressed 
in terms of integrals containing the crystal direct corre¬ 
lation function. 


v(q) = j <fVi Jd d r 2 n(r 1 )n(r 2 )t 


—*q-ri2 


x / jOng) 
V »( r i) 

v + 0(q 2 ), 
k B T 


- c(ri,r 2 ) 


Pot (q) = ^y~J d ^J ^V2c(ri,r 2 )(l-e- iq ri2 ) 


(14a) 
(14b) 

X _ e -*q r 12 N 


x n(ri)V a n(r 2 ) 
ip a pqi3 + 0(q 2 ), 

k B T 
V 


(14c) 

(14d) 


A a ^(q) = -jy- f d d r 1 J d d r 2 c{ ri,r 2 )(l-< 

(ri)) (v,gn(r 2 )) 

0(q 4 ). 


,-*q-ri2 


t 


[y a n 


(14e) 

(14f) 


The small wavevector limit and the index-symmetries 

Pal 3 Pi3a ^ad ^afi'yS ^Pa'yS X afi5~i X -ySa(3 are 

discussed in detail in Ref. |2- The explicit integrals are 
given in Eqs. (271, (301 and (32) below, where also crucial 
steps in their derivation are recalled. The connection of 
the elastic coefficients to thermodynamic derivatives will 
be established in Eqs. (39) and (431. 

The obtained set of equations (131 is best presented in 
matrix notation 


Vk B T5i-i= 


/ (5n* Sn) 

< 



(Su*Juf 3 ) J 


(n*(q) 

p*( q) A 

\pi 3(q) 

k 



kj 

(15) 


with Latin indices i = 0, a. The matrix of correlation 
functions of the macroscopic variables is thus given by 
the inverse of the matrix of the generalized constants of 
elasticity 


( < 'It > 

-i'S-Sup) 


(5u* a 8u 0 ) 


= Vk B T 


( Kq) 

p%{q) 

\Pa( q) 

Aa/3(q) 


-1 


(16) 
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The single matrix elements corresponding to the 
wavevector-dependent correlation functions ar^® 1 


when defects are neglected. Thus, it correctly captures 
the ideal crystal limit. Eq. (15) is transformed into 


A OlQ 




- W 

V 


{ ^) =V k B T ('- + <& 

TIq \U V L 

(17a) 

= Vk B T (u - _1 = Vk B TK~\ 

(17b) 

(5u*Ju p ) = Vk B T (\ af) - = Vk B THJ 

(17c) 


Vk B T (x~} + , (17d) 


,Sn* 
' no 


6up) = Vk B T 

= Vk B T (-JTV* A"*) 


~(S <^> = Vk B T 

= Vk B T (-A . 


(17e) 

(17f) 

(17g) 

(17h) 


The second line of each expression is a non-trivial alter¬ 
native, which is here given for completeness; it is based 
on the algebraic Woodbury identity. 

We thus reached our first goal of expressing the cor¬ 
relation functions of the coarse-grained variables, hydro- 
dynamic density and displacement vector field, in terms 
of integrals containing the functional derivative of the 
free energy with respect to density, viz. the direct cor¬ 
relation function. Let us note in passing that transla¬ 
tional symmetrjEI yields the expected g-divergences or 
g-dependences of the correlation functions in the limit 
q —> 0. In particular, (5u* a 8u,p) oc 1/g 2 follows from 
A a/3 (q) oc q 2 . 


2. Including defect density 

Although the relation between the constants of elas¬ 
ticity and the fluctuations of the coarse-grained fields is 
complete, it is worthwhile to consider a second set of vari¬ 
ables. So far the displacement field u a appeared in two 
different ways. It contributes to the coarse-grained den¬ 
sity, but it also appears as broken symmetry variable. In 
this section we introduce the point defect density c in lieu 
of the coarse-grained density, and keep the displacement 
field. 

The correlation functions of the coarse-grained den¬ 
sity and displacement field are easily transformed into 
correlations between the fluctuations of the point defect 
density and the displacement field using the definition 
Eq. §>. The set of variables {dc(q), <5u a (q)} may be 
more appropriate to describe an experiment when few 
defects are present and <5c(q, t) can be measured easily. 
It allows one to take the limit of vanishing defect den¬ 
sity and thus it is a natural set of variables to be used 


VksTSu = 


(<5c*<5c) 

«n 


( S <t 0 ) 

(■ Su* a 5u a ) 


ik 


v * (q) 

n 0 V s (q) 

n 0 V:(q) 

Ks( q) 


kj 


(18) 


The combination of the constants of elasticity appearing 
here is directly connected to the hydrodynamic equation 
of the momentum density expressed in terms of point 
defect density and displacement held as hydrodynamic 
variable^- 1 . There, the time derivative of the momentum 
density couples to the displacement held via the negative 
of 


A a ^(q) = A Q/3 (q) - iq a fj~p{q) + iti* a {q)q/3 + g Q ^(q)g/3- 

( 19 ) 

The coupling to the point defect density is given by the 
negative of 


K*(q) = — U* a ( q) - iq a v{ q)) • 

n 0 V I 


( 20 ) 


The individual matrix elements of the correlation func¬ 
tions in terms of n(q), Va(q), and A a/ g(q), and the limit 
q —> 0, may be determined according to the steps in 
the previous paragraphs. As the results can be obtained 
from Eqs. by straightforward replacements, identi¬ 
fied from comparing Eqs. (15) and (18), they will not be 
repeated here. 


B. Inverse relations 

Equations 0 predict the fluctuations of the macro¬ 
scopic coarse-grained density and displacement held 
based on the generalized constants of elasticity obtained 
from the direct correlation function and thus the free 
energy. Experimentally, the inverse relations are of in¬ 
terest: expressing the generalized constants of elasticity 
of the crystal in terms of measurable correlation func¬ 
tions. Two different sets of correlation functions can 
be obtained from experiments. Either displacement held 
and coarse-grained density fluctuations can be recorded, 
or displacement held and defect density. For reference, 
we provide the inverse relations for both cases in this 
section. 


1. Including coarse-grained density 


The matrix equation (15) can be inverted in order to 


find the elastic functions i'(q), A i a (q), and \ a p(q) in 
terms of measurable fluctuation functions. The inverse 
relations read: 
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Kg) 

Vk B T 


Ag^(q) 


Ma(q) 


,671*671 


n\ 


-1 


0 


'T.g rig 


ns 


n 0 


n 0 


/ dn^ d77_ 1 <^^—)) =fl \ 

V ng n 0 H n 0 ) 

(< 6 < Su fi ) - ( K - K ^)-^-^))" 1 = 5 - 

V no nn n n / p 

Sus)(5u* s 5up)~ 1 , 


A77 Aon 

{5u* a 5up)~ l + (Su* Su J )~ 1 (Su* — )R~ 1 ' 


n 0 


n 0 


S-t(Su}-)( 


Sn,,5n*5n 


n 0 
-1 / 


«n 


-1 


<5n. 


(Su* a Sup) (Sul — )i? 

no 


" 1 . r * <5n Sn*5n 
{ Su } — ){ — 2 - 
p n 0 ni 


)- 1 (21a) 
(21b) 
(21c) 
(21d) 
(21e) 
(21f) 


We thus reached our second goal to derive relations which 
determine the generalized elasticity constants A a p(q), 
p a (q), and n(q) from measurements of correlation func¬ 
tions. 


2. Including defect density 


Replacing the total density with the defect density us¬ 
ing Eq. the generalized constants of elasticity can 
be connected to fluctuation functions which can be mea¬ 
sured at constant (possibly vanishing) defect density. 
The comparison of the matrices in Eqs. (151 and (18) in¬ 
dicates the straightforward replacements in Ecp (21). The 
dynamical matrix A a/ g(q) determines the wave equation 
of the momentum density, and its eigenvalues give the 
(acoustic) phonon dispersion relations. Because the re¬ 
sults follow from straightforward replacements, they will 
not be given explicitly here. 


IV. FREE ENERGY AND THERMODYNAMIC 
RELATIONS 


In order to obtain the thermodynamics derivatives, a 
consideration of the free energy appears useful in cases 
where the connection to the small wavevector limit of 
the correlation functions is not established or under 
debat^S^H In this section, we will coarse-grain the free 
energy functional of density functional theory in order to 
derive the thermodynamic derivatives. This purely equi¬ 
librium statistical mechanics approach supplements the 
dynamical one based on projection operator formalism 
in Ref. [2J. Importantly, the wavevector dependent cor¬ 
relation functions of the coarse-grained fields of elasticity 
theory and the thermodynamic elastic free energy of real 
(viz. defect containing) crystals are then obtained from 
a single microscopic starting point. 


A. Coarse-grained free energy functional with 
elastic fields 

The second order change in free energy AT due to a 
deviation Sp( r) in the average density distributi on fro m 
the periodic crystalline equilibrium density n( r) iJSSlST 33 

AJC = LbT_ JI - c(ri,r 2 ))5p(ri)£/9(r 2 ), 

( 22 ) 

where c(ri, r 2 ) is the direct correlation function of a peri¬ 
odic crystal. Note that this quadratic functional contains 
the direct correlation function as single input and thus 
the identical information as used in the correlation func¬ 
tions approach of the previous Sect. [Tl} 


1. Including coarse-grained density 


We start from the representation of the microscopic 
density fluctuation in terms of displacement field and 
coarse-grained density, Eq. ([ 5 ]). We assume that an anal¬ 
ogous equation holds also for the averaged (macroscopic) 
densities. In this way we get a change of the average den¬ 
sity due to non-vanishing displacement field and average 
coarse-grained density, 

Sp(r) = —<5u(r) ■ Vn(r) + n(r ) ^ n ^ , (23) 

n 0 

We shall emphasize that while Sp( r) varies on the spatial 
scale of the crystalline lattice, the coarse-grained density 
varies far more smoothly and contains wavevector con¬ 
tributions only from the first Brillouin zone: 


Sn( r) = f 

Jv 


1 st BZ 


d d q 

(2n) d 


6n(q) . 


Using Eq. (23) we obtain the following expression for 
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the product of density changes 
6p(r 1 )Sp(r 2 ) =5u a (r 1 )6up(r 2 )\' a n(r 1 )V pn(r 2 ) 

" -V-' 

(!■) 

, n(r 1 )n(r 2 )5n(r 1 )Sn(r 2 ) n(r 2 )<5n(r 2 ) 

H- 2 -ow a (iq)V a n(iq)- 


Our subsequent calculation is analogous to that of 
Masters^ and is equivalent to the discussion of the sur¬ 
face tension in 135] . We will in detail describe the calcu¬ 


lation originating from the first part of Eq. (24), which 


leads to the elastic tensor A, and then summarize calcu¬ 
lations originating from the other parts. 


n 0 


(2.) 


(3.) 


n(iq)<5n(iq) 


n 0 


<5u a (r 2 )V Q n(r 2 ). 


(24) 


(4.) 


In the following, we substitute the four parts of Eq. (24) 


into Eq. ( |22| ). We denote the resulting expressions AJq.j ), 
where i = 1,..., 4. We then re-write these expressions 
using the equation 


In the expression for .) one expands Sup(r 2 ) 

around rq, which is valid for a short range (in ri 2 ) di¬ 
rect correlation function. The zero order term vanishes, 


because of (25) and the first order term does not con¬ 


Vqn(r) 

n(r) 


= J dVc(r, r')V Q n(r'). 


(25) 


tribute due to the symmetry c(ri,r 2 ) = c(r 2 ,ri). Since 
the hydrodynamic variable 5u(r) is slowly varying, one 
obtains an expression which is quadratic in V(5u(r) as 
leading contribution 


A^(i.) = 


- c(ri, r 2 ) b5M a (ri)i5M^(r 2 )V Q n(ri)V / 3n(r 2 ) 


JJ d d, ri d d r-2 ( 

JJ d d r 1 d d r 2 Vo,n(r 1 )c(r 1 ,r 2 )Vpn(r 2 )5u a (r 1 )(5up(r 1 )-5u0(r 1 ) + y i Sup(r 1 )r 1 2,^ -l i V'y\7sfiup(r 1 )r 1 2,- / r 12 ,s 


S(ri2) 

n( ri) 


= -//d c V 1 d'V 2 


k B T 


V Q n(ri)c(ri,r 2 )V^n(r 2 )ri2, 7 ri2,,5 


=0 symmetry 
V 7 i5it Q (ri)V55u / 3(ri) 


d^'Xap^tV-fSua (r)\7$5up (r), 


k T f f 

XapjS = / / d £ Vid t V2V Q n(ri)c(ri,r 2 )V /3 n(r 2 )ri2 l7 ri2,5 . 


(27) 


In the last line of Eq. (26) the separation of spatial scales With 


was used in order to replace the slowly varying local elas¬ 
tic coefficient given by the contents of the square bracket 
on the third line of Eq. ( [26] ) by the macroscopic constant 
of elasticity X a p-yS- We emphasize that the expression 
(27) agrees with the one obtained in the framework of 
hydrodynamic equations of motion^, which was repro¬ 


duced in Eq. (14 


For the second term of the free energy, Sn(r 2 ) is ex¬ 
panded around iq and, as hydrodynamic variable, as¬ 
sumed to be slowly varying 




x [n(r 1 )<5(r 12 )-n(r 1 )c(r 1 ,r 2 )n(r 2 )], 

(28) 


/M^)’ 


(29) 


k B T 

V 


J J crtqcrtq(n(iq)< 5 (ri 2 )— n(iq)c(iq, r 2 )n(r 2 )) . 

(30) 

The third and fourth part yield with the same arguments 

AJ (3i+ 4 .)=- [ d d r p a p X7pSu a (r), (31) 

J n o 

k T f f 

= / «Hq / d < V 2 n(ri)V a n(r 2 )ri2,/3c(ri,r 2 ) 

(32) 

Summarizing, we obtain the following expression for 
the free energy change 


A T = 


j d\ + Ca W v(r)«ii(r) 

J ' n 0 / 

f ,d M r ) 

- / d\p a p - u a , ?(r) 

J no 


(33) 



































Expression (33) involves the symmetrized linear strain 
tensor u a p(rj= a dup(r) + \7 p5u a (r)) and the Voigt- 
symmetric elastic constants C™g s = X a -yp$ + Xp ia $ — 
A a pjs- Both combinations reflect the rotational sym¬ 
metry as only symmetric combinations of strain enter 
into the elastic energy and the tensor of elastic con¬ 
stants obeys a number of symmetry relations. Their 
prooP is based upon the rotational analog of the LMBW 
equatioriSS 


ri x V (1) lnn(ri) = f d d r 2 c(r 1 ,r 2 ) (r 2 x V^n(r 2 )j. 


(34) 


This free energy functional is a superposition of inde¬ 
pendent terms each containing the square of the Fourier 
transformed coarse-grained fields at one specific wavevec- 
tor. Often one connects such quadratic free energy 
functionals with a probability distribution for fluctua¬ 
tions of the coarse-grained fieldspl P[<5n(q), <$u(q)] oc 
exp {—AP/fegT}. In the present case, this would yield 
the wavevector-dependent correlation functions © as 
statement of the equipartition theorem resulting from 
this Gaussian approximation. 


B. The thermodynamic elastic free energy 


We thus arrived at our third goal, to derive the general 
elastic free energy functional of real crystals containing 
the coarse-grained macroscopic fields. Let us add that 
the above expression for the free energy also determines 
the constant Co = 0 in Eqs. (89), (90), and (92) of 
Ref. [2], which could not be determined from the hydro- 
dynamic equations considered there. 

2. Including defect density 

The free energy in terms of the defect density <5c(r) 
and the displacement field <Ju(r) is obtained from Fourier 
transforming ansatz ff and Eq. © into real space: 

dp(r,t) = —V ■ [n(r)Au(r, f)l — ^<5c(r, t). (35) 

n 0 

Following the steps of the previous section one arrives at 
the coarse-grained free energy including the defect den¬ 
sity: 

AT = + 2(vS a p + p a pj ^^-u a p(r) 

T a.P'yS ~^~^d a pSryS T Papd^yS T d a pP'yS^ ttap (y)'U’jS (r ). 

(36) 


The thermodynamic free energy corresponds to homo¬ 
geneous fluctuations, viz. the coarse-grained fields eval¬ 
uated at q = 0. It can handily be obtained from the 
explicit free energy functional in Eq. (|33|). The result 
shall be given using the Voigt notatiorPT^m three dimen¬ 
sions), because this form appears convenient for explicit 
model calculations later on. Quantities in Voigt notation 
carry lower Latin indices 1 < * < 6, where u, denotes 
the independent elements of the symmetric strain ten¬ 
sors u a p. For 1 < i < 3 the relation iq = u a _p holds 
with (a,/3) = {(1,1); (2, 2); (3,3)}, while for 4 < i < 6, 
Ui = 2u a: p holds with (a, 0) = {(2, 3); (1, 3); (1,2)}. For 
spatially constant fluctuations (to be indicated by sub¬ 
script q = 0 where otherwise unclear), one obtains in 
obvious notation as a quadratic form: 



The thermodynamic free energy is a quadratic form given 
by a 7 x 7-matrix of elastic coefficients, where in Voigt 
notation the elastic matrix is C l3 = C a p 1 $ for 1 < i,j < 6 
with the index correspondences given above. 

1. Connection to elastic coefficients and variances 


This gives the relation between the elastic coefficients at 
given defect density C c in terms of the corresponding 
coefficients at given total density, C n , namely: C^g s = 
^aP'y5 d" ^dapdryS T Papd~yfi 4“ dapp~yS • 


3. Gaussian probability distribution function 


The harmonic free energy Eq. (|33j) can be written in 
a more compact form with the help of the 4 x 4-matrix 
of elastic coefficients introduced in Eq. (15). Fourier- 
transformation leads to 


^=1 


Sn*( q) 
, no 


d d q 

(27r) d 

, <5n*(q 


v -ip-ysqs 
ip a pqp C™ MS qpq& 


(37) 


<fa(q) 

n 0 

6u 7 (q) 


Thermodynamic derivatives can now easily be evalu¬ 
ated and lead to the parameters already introduced in 
Eq. (14). They follow from the Gibbs fundamental form 
of the free energy density / = F/V « AF/V, where the 
quadratic expression (38) suffices in order to obtain the 


second order derivatives of interest. 


aV 

dn 2 

d 2 f 

dndu a p 

d 2 f 


dp 

dn 

dp 


du a pdu lS n durys 


du a p 

dh a p 


= i//r 


dh, 


Vi 


ot/3 


dn 


(39a) 

= -ptap/no, (39b) 


— (-'aP'yd — ^a'ypS V A P'yaS X a p^5- 

(39c) 


These relations identify the elastic parameters of our 
approach as thermodynamic derivatives. They already 
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use the familiar intensive variables, chemical potential /i 
and str ess t ensor h a p in order to familiarize with later 
relations 1 3 These variables will be introduced and dis¬ 
cussed in Sect. IIV Cl below. Let us note that these cal¬ 
culations supplement the derivation of the thermody¬ 
namic relations in Ref. [2] (recalled in Eq. ©) , where 
the equivalence of the hydrodynamic equations was used. 
The thermodynamic free energy thus takes the form: 


AT = — ( — 

2 V n° 


Ui 


_ dhj dhj 

n 0 dn dm 


Sn \ 

™ j (40) 


Where, in Voigt notation the stresses correspond to h, = 
h a p for 1 < i < 6. 

This compact expression is a convenient starting point 
for evaluating the thermodynamic covariances and sus¬ 
ceptibilities which enter elasticity theory. The isothermal 
compressibility and the defect density susceptibility will 
be obtained in the next Sect. |IV C[ In order to prepare 
for this, first the second moments of the fluctuations of 
the thermodynamic variables shall be obtained. These 
are connected to the thermodynamic derivatives using 
the thermodynamic formalism. Because the inverse of 
the Jacobian matrix is equal to the Jacobian matrix of 
the inverse function one obtains 


v 

-Hi 


~Hj N 

/-in 

y 

r-( 

,-»« s 

du 

~ n °duf 

dhj 

duj 


-( 

^ 1 dn 

n q d/i 

_ 1 dn 

\ no dhi 

1 duj 
no d/i 

diLj 

dhi 



1 

Vk B T 


,Sn v 
n j 'n 0 3/ 

(' u i W> ( u i u j) 


(41) 


< 2=0 


In the last step the fluctuation-dissipation-theorem is 
usecPS The variance of the total coarse-grained den¬ 
sity variation is thus obtained from a simple matrix 
inversioii 30 : 


f 6n5n 
\ „2 / 


9=0 


= Vk B T - + 


1 , Hi 


v v L 


/—in 


[liflj 


V 


(42) 


= Vk B T(v-Hi(C? j )~ 1 H j 
= Vk B T[v~H a p 


H~yS 


2. Including defect density 


In a similar manner an expression for the defect den¬ 
sity fluctuation can be obtained. Starting from the free 
energy functional in Eq. ( |36| ) and considering homoge¬ 
neous variations (viz. at q = 0), one recognizes that the 
relevant thermodynamic derivatives are now given by 


av 

dc 2 u a 

d 2 f 

dcdu,, 


dn 

"dc 

dn 


— v/r 


'O’ 


(43a) 


a/3 


du, 


af3 


da, 


aj3 


dc 


d 2 f 


duafidu-fti 


= (v8 a p + Hap) /n 0 = Hap/ n 0’ (43b) 

— ^a.p'yd — afi'y8 "b Hctft&'yS V SafiH'yS 4“ , 

(43c) 


where the stress tensor a a p was introduced, which will be 
discussed in in Sect. IIV Cl below. Also the abbreviation 
fi c was introduced. Thus, using the fluctuation dissipa¬ 
tion theorem again, the matrix of total thermodynamic 
variations is given by 


l i(te Uj ) 

n£ 'n 0 ■?' 


= Vk B T 

< 2=0 


= Vk B T 


— 1 dc 

7T.q 8/1 

1 dc 
no d&i 

v Hj 
Hi C9. 


\ 

n 0 du | 
duj I 

d^i ) 

) (44) 


C/I and Hi are the tensors from Eq. (43c) and Eq. (43b) 
in Voigt notation. This leads to the correlation of the 
defect density fluctuations 


, ScSc 
'Vk B Tn 2 


< 2=0 


(; 




— (v ~ HapiCapjt) V75) 


H c ' 

V V L 
-1 

■o) 

-1 


/-in 


llitlj 


(45) 


As in Eq. (42), the second line followed from a Woodbury 


identity, and the usual notation is used instead of the 
Voigt one in the last line. This covariance of the number 
of point defects, will be connected to a compressibility¬ 
like expression k c below. 


C. The isothermal compressibility of crystals 


where the second line follows from a Woodbury identity, 
and the usual notation is used instead of the Voigt one in 
the last line; see Wallac^l and the Appendix [A] for the 
proper interpretation of the inverse fourth-rank tensor. 

We thus derived the second moment of the particle 
number fluctuations from DFT. We started from the 
same free energy functional as was used in the deriva¬ 
tion of the wavevector-dependent correlation functions 
summarized in Eq. Thus, in Sect, [v} both results 

can be compared in the long-wavelength limit. 


In this section the general expression for the compress¬ 
ibility of a real crystal is derived from a thermodynamic 
consideration, details are given in Appendix [A] The sit¬ 
uation described is one in which no pre-stress is applied 
to the crystal in equilibrium. 

The definition of the isothermal compressibility of a 
fluid reads 


1 dV 

yjj) 


N 


( 46 ) 
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where p is the pressure. For a crystal, the question arises 
how this has to be generalized to describe the additional 
degrees of freedom. The infinitesimal change of the free 
energy of a crystal at constant temperature (dT = 0) 

dF = —pdV + pdN + h a pdU a p, (47) 

includes a term with a stress tensor h a p at constant 
volume V and particle number N times an extensive 
strain tensor U a p = Vu a p. The work done is SW = 
f h a pSu a pdV with the symmetrized linear strain ten¬ 
sor u a p = |(V a U /3 + Vpu a ). The chemical potential 
is denoted by p, and the particle density will be denoted 
n. While this ’first law of thermodynamics’ for a crys¬ 
tal is familiar from textbookJO the coupling of strain 
and density fluctuations complicates the interpretation of 
the stress tensor h a) 3 , which calls for a discussion before 
addressing the compressibility. Taking a ca nonica l N- 
particle system and straining it infinitesimally 37 * 3S I leads 
to the stress tensor t a p = TAG- at fixed T and N. (It 
can be also obtained from averaging the Irving-Kirkwood 
microscopic stress tensor.) Because the volume V varies 
in this procedure, the two stress tensors differ by a scalar 
term £1211 — (p— d a p- The compress¬ 

ibility for a periodic crystal shall be understood as the 
derivative at constant h -stress tensor, because it then 
measures the change in particle density with chemical 
potential, 


This expression for the isothermal compressibility of a 
general crystal generalizes results obtained for high sym¬ 
metry crystals 1 " Hence, together with Eqs. (50) and 
(42) and Sect. Ill A we achieved our main goal to es¬ 
tablish the general connection between the isothermal 
compressibility of non-ideal crystals and the correlation 
functions of the fields of elasticity theory. The connec¬ 
tion is derived from microscopic DFT. See Appendix [A] 
for an alternative formulation of Eq. (51) derived within 
the thermodynamic formalism, and corresponding to the 
first line of Eq. ( |42| . (Eq. (51) corresponds to the second 
line of Eq. (42).) 

If the coupling between strain and density fluctuations 
in the result for the isothermal compressibility in Eq. (50) 
is neglected, the second term vanishes and the compress¬ 
ibility k is given by kT 1 = v = , which plays 

n «c«3 

the role of the inverse bulk modulus at constant strain. 
While in regular solids, the coupling between strain and 

• r 

density in the free energy, p = dn g u , cannot be neglected 
and this approximation fails, see Sect. VI for a system 
where it holds well. In order to dissect the contributions 
to the compressibility in detail for more regular crystals, 
transforming to defect density is required. 


2. Including defect density 


1 dV 

Vhp 


N,h n 


1 dn 

rig dp 


(48) 


where we used Maxwell and Gibbs-Duhem relations de¬ 
scribed in Appendix [Aj They lead to the Gibbs funda¬ 
mental form of the free energy density / = F/V which 
was already anticipated in Eqs. (39), namely: 


df = pdn + h a pdu a p . 


(49) 


Also, the calculations for determining k have already 
been done. Equations. (41) and (42) immediately give 


the isothermal compressibi 
density fluctuations: 

1 


ity as variance of the total 


,5n5n, 


Vk B T x nl 


9=0 


1. Including density 


While the result for k in terms of the elastic coefficients 


(viz. Eqs. (42) and (50)) is useful for explicit evaluations, 


and will be used in Sect. [VI] below, a relation connecting 
it to thermodynamic derivatives is desirable and would 
take the form expected in the thermodynamic formalism. 


Using the relations (39) in order to replace the elastic 


coefficients in Eq. (|42|)7we find 

'dp 
>.dn 


1 


K = 


dh, 


a/3 


dn 


/ dh^s \ 

V du a p nJ 


If one considers the set of independent variables with 
the defect density c instead of the coarse-grained density 
n with Eq. (ro|) simplifying to 


dn = —nodu aa — dc . 


(52) 


the manipulations leading from Eq. (1501) to Eq. (511 have 


to be adapted. The compressibility is given now in terms 
of derivatives at constant defect densitj®. The stress ten¬ 
sor (7 a p (with cr a p = h a p — nopS a @) and the chemical 
potential p now are functions of the strain tensor and 
the defect density combining to the Gibbs fundamental 
form of the free energy density - 1 df = —pdc+cr a pdu a p. 
The relevant thermodynamic derivatives are given by 
Eqs. (43) (see Eq. (36) for the free energy density), which 


(50) need to be used in order to replace the elastic coefficients 


in Eq. (42). This leads to: 
' 3 . dn 


n = — \ n, 



-1 


+ \5 a p- ( 


-1 


dp 


du. 


a/3 


dcr-fS 


(dp 

y 1 dp 

c dc 

U a /3 

{dc 

u a ,3 J du a p 





2 dp 

1 dh^s 

y 1 

- - t 

dn 

U a pJ 

(51) 

X 

' da a p 

du^s 



-1 N 


c , d(T a p 

Oafi + - 

/3 OC 


( 53 ) 


- n\ 
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dcr at 3 


-no ^ 7 5 no ~ 

OC U~yS OC 


X 


f 2 d d 

{ n °Tc 


U~y§ 


X I ^ d&ryS 

0-y5 H - no 


U~yS 


U a /3; 


(54) 


An interesting limit is now the vanishing of the cou¬ 
pling between the defect density and the strain field, 
g e g M J = 0; see Eq. (|43bj). This yields two indepen¬ 


dent contributions to the compressibility which shall be 
denoted n° in this approximation 


= i," 1 


+ (CaB~/s) — V + ^ (Qj) 


E 

*,i=i 


\-i 


(55) 


The first contribution v _1 is due to the fluctuations of the 
defect density, whereas the second one 6 )~ x 5 a ^8ys 

is due to independent fluctuations of the strain tensor, 
which agrees with the known result for a perfect crystal 
without external strairP^. 


D. The isothermal defect density susceptibility 

Varying the chemical potential changes not only the 
average density but also the defect density. The deriva¬ 
tive of the defect density with respect to ji can be ob¬ 
tained analogously to Eq. (50), and a thermodynamic 


susceptibility akin to the compressibility can be defined 


— 1 dc 
7T-0 dfx 


. Sc5c . 
Vk B Tnl 


9=0 


(56) 


The explicit result for k c in terms of the elastic coeffi¬ 
cients is given in Eq. (45) and in terms of the derivatives 


from Eq. (43) is given here: 


k = 


u 


dfj, 

dc 


da. 


aft 


u a 0 


dc 


/ da-ys_ \ 
\du a p cJ 


u a p 


days 


dc 


u a fj) 

(57) 


Connecting the isothermal defect density susceptibility 
to derivatives of the density appears useful in order to 
obtain it e.g. from computer simulations. Starting from 
the definition of k c in Eq. (56), the Eq. (52) leads to 


-1 dc 
o d ! 1 


n\ 


— 1 dc dn 


Hq dn dfi 


— 1 dc du, 


n, 


1 dn 


+ 


1 du„ 


at/3 


du a/3 dfi 


no <9/i 


(58) 


The first term on the right hand side is a thermody¬ 
namic susceptibility at constant er-stress tensor, which 
bears similarity to the isothermal compressibility: 


n\ 


1 dn 
o 


(59) 


Yet, Appendix[B]will show that this specific susceptibility 
vanishes in the limit of an ideal crystal, and thus does not 
play the role of a compressibility in solids. The second 
term can be reformulated using Eq. (44), and the result 
can be rearranged to give: 


K a = K c 1 - 


days 

dc 


da n 


l V dUyS 


( 


n 


(60) 


= k c (i - ^(c^j- 1 ) = * c (i - ^(c&r 1 )* 

i =i 

where in the last equality the thermodynamic derivatives 
from Eq. (431 were used, and the result transferred in 
Voigt notation. The difference between k c and n a , which 
both are derivatives at constant er-stress tensor, vanishes 
in cases where strain and defect density fluctuations do 
not couple (viz. /z c = 0). In the general case, density and 
(the negative of the) defect density vary differently with 
chemical potential at fixed a. 


V. SMALL WAVEVECTOR LIMIT OF THE 
STRUCTURAL FUNCTIONS 

So far we considered correlation functions and the 
isothermal compressibility of crystals. In this section we 
bridge the gap between the density correlation function 
and the compressibility, and point out the subtle differ¬ 
ence between the two expressions. In the second part of 
this chapter the so called generalized structure factor is 
discussed. 

In order to understand the connection to the compress¬ 
ibility, the g-dependence in the limit q —> 0 of the corre¬ 


lation function of the coarse-grained density (17a) needs 
to be discussed in detail 


Sn*Sn 


)=u (q) + v (q)/x*(q) 

-1 


Vk B Tnl' 

x (-Wq) - Ma(q)^ _1 (q)/4(q)) M 7 (q)^ _1 (q) 

fJ"y5Q5 


9^0 1 ^ flgpqp 
V V 


(A, 




0:7 eq. 


)<itq<t> 


(61) 


Here we used the known small-wave vector expansions of 
the elastic coefficients, which were defined in Eqs. (14). 


They follow from DFT relations expressing translational 
and rotational symmetry®. Noting that only the sym¬ 
metrized combinations in a o 7 and e f> of the term 
in square brackets are relevant, and with the help of 


Eqs. (39) this expression can be further simplified to 


Sn*Sn 9 -h) 1 ^a/ 39/3 


v Vk B TnV 


(r ,n 

V'~'a:e 70 )Q.eQ.(j) 


-1 


M 'y8Q8 
(62) 


This expression would agree with the thermodynamic 
one (50), if the factors of qpqs canceled g c g^. That the 
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limit q —> 0 is not that simple can be seen even for highly 
symmetric crystals. For a cubic crystal, the correlation 
function yields different results in the small q limit (62) 
depending on the direction of q relative to the unit cell. 
And for the hypothetical model of an isotropic crystal, 
the small q limit (62) is direction independent, but dif¬ 
fers from the thermodynamic value from (501. The latter 
simplified case, allows to identify the origin of the discrep¬ 
ancy and will be studied in detail in the next section. 


A. Perfect crystal embedded in a matrix 


To study the difference in more detail, it is, as a first 
simplification, more convenient to look at the simpler 
problem of a perfect crystal. In this section we also use 
the more familiar expressions of elasticity theory. The 
connection to the terms used so far is given in Appendix 
|B| For a perfect crystal the correlations of the displace¬ 
ment field is given by the (inverse) of the dynamical ma¬ 
trix D a p( q) (for particles with mass m) 


{5u* 5up) = —— D p(q). 

mriQ p 


(63) 


The coarse-grained density fluctuation for a perfect crys¬ 
tal is (5n(q, t) = —inoq a Su a (q,t) and the dynamical ma¬ 
trix is related with the elastic constants 1 ^ via D a7 (q) = 
C'ap-ysqpqs- Thus for the coarse-grained density correla¬ 
tion function we get 

Sti* 5n 

( Vk B Tn 2 ) = d.a.D a p i ( l)qi3/{ mn o) = Qa{C' ae / 3 ( fiq e q ( f l ) qp 

(64) 

which shows the same problem in the limit q —» 0 
as arises in Eq. (62), when compared to the thermo¬ 
dynamic compressibility of an ideal crystal K lc = 
{C~p ri& )8 a p8 1 s (contraction of the inverse of the matrix 
of elastic constants). For an isotropic crystal the elastic 
tensor simplifies to the two Lame coefficients C a pjs = 
X5 a pSjs + p(S ai dps + S a sSf 3 7 ). Thus, the compressibil¬ 
ity is (/c lc ) _1 = A + whereas the correlation function 
yields X+2p (which corresponds to the longitudinal speed 
of sound). 

To show the origin of this difference we consider an 
isotropic (ideal) solid for which the so called fundamen¬ 
tal solution of elasticity is known. Other symmetries with 
known solutions are hexagonal and pentagonal. The 
corresponding problem in two dimensions can be found 
in [42j- We consider a three dimensional sphere with 
volume V B embedded in a spherical matrix V of the 
same isotropic material. The radius Rb of the embedded 
sphere is increased R B -A R B + Ar and the surrounding 
matrix is compressed. To determine the displacement 
field and the elastic energy of such a deformation one 
has to solve the equation of elastostatic theory, which is 
the vanishing of the divergence of the stress tensor, or in 
terms of displacement field 


^ — 0 - 


(65) 


The solution is a sphere with increased volume V B +AV B . 
The only non-vanishing displacement field is (homoge¬ 
neous dilatation) 


Su r 


Ar ife r<R B 

Ar(Aa) 2 r>R B 


This yields for the total elastic energy 


E= v B( Av B y 
2 V V B J 


+ g A 1 ) 



VbY 
~V ).' 


( 66 ) 


(67) 


The first part is due to the stretched sphere and the sec¬ 
ond contribution is from the surrounding matrix. Thus, 
depending on the ratio V y the relevant combination of 
elastic constants changes from A + | p for (^f- -A 1) to 
A + 2/i for —> 0). In the limit of vanishing shear 

modulus p. the difference vanishes. Thus, for a fluid it 
doesn’t matter if one determines the volume fluctuations 
of a small sphere in surrounding fluid (of the same kind) 
or if one looks at the global fluctuations of the whole 
system. 

It is worthwhile to note that the same ratio between 
these two combinations of Lame coefficients appears in 
a related context. In Eshelby’s study^ of an inclusion 
in a matrix of elastic material, the so called constrained 
strain is given by the stress-free strain u ^ 


u 


C 

aa 


A + 

A -(- 2p 


( 68 ) 


This calculation has recently been extended to atomistic 
inclusionJSI, which could be used to test approximations 
in the present DFT approach. An a nalogou s problem is 
a polar fluid in a dielectric mediunP^ESEH. There, the 
susceptibilities show a directional dependence due to the 
dipolar interaction, and a different combination of dielec¬ 
tric constants is relevant depending on the surrounding 
medium. 


B. Generalized Structure Factor 

There is a further aspect which differs the relation be¬ 
tween the compressibility and the correlation of the den¬ 
sity fluctuations of a fluid and a crystal. There is a dif¬ 
ference if one looks at the elements of the generalized 
structure factor which contribute to the compressibility, 
i.e. whether those are different from S’ g _ 0 (q -A 0). 

We recall that the generalized structure factor is de¬ 
fined bjP 9:: 

S g ( k) = i J (fVr J d% (5p(r 1 )<5p(r 2 )) e- ,g R e^ k Ar 

= g/2 + k)Sp{g/2 - k)) , (69) 

(with R = (ir+r 2 )/2 and Ar = iq — r 2 ) and its 5p(g + 
q) element is measured in a scattering experiment 1 36 . 
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In a liquid, where translational invariance dictates 
that only Sp(k ) is non-vanishing and isotropic, its 
connectioi i — I ^ to the compressibility is well known 
Sq (q —>• 0) —> ri^ksTn. To convince oneself that such 
a connection does not hold in a crystal, the definition of 
the coarse-grained density Eq. ([8]) can be used to derive 


Sn*Sn 


= jp 5Z ? v( ( V‘(g' + q)Mg + q)K 

0 g,g' 

= (2*) fl 
" M 3 


f - s- q)£(g - - §> 


g,g' g 

U g,g' 


g-g 


g + g' 


q n g , 


(70) 


where the left hand side becomes k for q to 
zero in the fluid case. Clearly, every element of 
(Sp*( g' + q)6p(g -f q)) is involved, not just the one with 
vanishing reciprocal lattice vector g = g' = 0. Also 
the correlation of coarse-grained density fluctuations is 
given by a combination of generalized structure factors 
S s - g '(—{g + g')/2 — q) in the limit q —> 0 and not just 
by Sg =0 (q —> 0) as for a fluid. Although the possibility 
that the right-hand side of the last equation is indeed 
the compressibility cannot be ruled out, it seems rather 
unlikely. 


VI. AN EXAMPLE: CLUSTER CRYSTALS 

To test the theory presented in the preceding sections, 
single component crystals of Bravais symmetry formed 
by spherical particles provide the closest cases. Large 
densities of local defects are desirable since the strength 
of the theory is its ability to account for the coupling 
of strain and defects densities. Additionally, a good ap¬ 
proximate DFT functional should be available. Recently, 
cluster crystals made from soft particles were discovered 
which satisfy these criteria and are thus ideally suited for 
testing the theory. 


A. Model and approximate density functional 
theory 

We consider a system of spherically symmetric par¬ 
ticles interacting via a purely repulsive, bounded pair 
potential. Following earlier studied®, we use a gener¬ 
alized exponential model of exponent four (GEM-4), 

$(r) = ee" ( - )4 . (71) 

The GEM-4-system shows several interesting properties. 
The finite upper bound of the potential allows cluster 
formation, i.e. the occupation of one lattice site by sev¬ 
eral particles. The GEM-4-system crystallizes in the fee 


and bcc phases, and at low temperatures it undergoes 
isostructural phase transitions between fee phases with 
integer occupation numbers per lattice site. At higher 
temperatures hopping of the particles between the lat¬ 
tice sites is possible and leads to a continuous, average 
occupation number. For the average density distribution 
of the cluster crystal the following ansatz is choserf® 

P(r)= «.0 i E^ R) ' < 72 > 

R 

with the occupation number n c , the inverse width of the 
(Gaussian) density distribution around a single lattice 
site a and the lattice vectors R. With this ansatz and 
an appropriate free energy functional one can get the pa¬ 
rameters n c , and a for given temperatures, and average 
densities through minimization of the functional. Then, 
the equilibrium state can be found by a direct compari¬ 
son of the free energies of each phase. For the description 
of the phase-diagram of the GEM-4, Mladek and cowork¬ 
ers found that a liquid-like mean-field approximation is 
appropriate® which leads to the simple expression for 
the direct correlation function c(ri,r 2 ) 

c(ri, r 2 ) = c(r) = —/34>(r) , with r = |r 2 — ri|. (73) 
This results in the following free energy functional 

F[p)=F id [p]+F ex [p], (74) 

F id [p\ = ^ j d 3 r[p(r) ln[p(r)A 3 ] — p(r)], 

F e*[p\ = \ J d 3 r 1 p(r 1 ) J d 3 r 2 $(ri,r 2 )p(r 2 ). 


Here A denotes the thermal de Broglie wavelength. By 
subtracting the free energy of the fluid from the crystal 
one, the parameter A can be eliminated without changing 
the position of the minimum of the crystal free energy 
functional. Similarly, dividing by the average density n o 
does not change the free energy functional minimum, but 
leads to a convenient expression 



(75) 


with the Fourier transformed potential Ju. As we alluded 
to earlier, the free energy functional (75) is to be mini¬ 
mized with respect to n c and a. The resulting, normal¬ 
ized free energy only depends on the single (dimension¬ 
less) thermodynamic parameter -. Thus, the fluid- 
bcc and the fcc-bcc phase transitions lie on straight lines 
drawn from the origin of the T — u 0 phase diagram. It 
should be noted that the free energy functional (751 is 
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minimized by the ratio instead of n r itself. Numer- 
ical minimization shows that varies only by about 
±3% in the whole solid phase, i.e. the system changes its 
density mainly due to changes in the occupation number 
and not due to changes in the lattice constant 


B. Compressibility and occupation number 
covariance 



FIG. 1: Compressibilities of the GEM-4 system in units of 
[noecr 3 ] -1 versus the reduced thermodynamic variable • 

While k is taken at fixed stress h a p, k c is taken at fixed stress 
a a p, and k° is the approximation neglecting the strain-density 
coupling introduced in Eq. (551. The approximation re « 1/V 
to neglect the strain-defect density coupling holds within the 
line thickness; see Fig[4] 


After minimizing the free energy and obtaining the av¬ 
erage density profile, the elastic coefficients from Eq. [l4j 
which are relevant for the compressibilities, can be cal¬ 
culated by straightforward integrations in the reciprocal 
space. The thermodynamic derivatives then follow from 
the relations in Sect. |IV C| Figure |T] shows three com¬ 
pressibility like quantities in all stable phases obtained 
from the mean-field DFT functional ([751. The compress¬ 
ibility k is taken at fixed stress h a p and describes the den¬ 
sity change with chemical potential /i. The susceptibility 
k c is taken at fixed stress a a p and captures the defect 
density change with /i. The quantities re° and \/v are 
approximations neglecting the strain-density and strain- 
defect density coupling, respectively. In reduced units, 
the thermodynamic derivatives change little throughout 
the complete stable fee phase. In the bcc crystal, de¬ 
fect fluctuations grow appreciably with increasing tem¬ 
perature. The fluid is less compressible than the solids, 
as follows from the familiar expression of the isother¬ 
mal compressibilitj^H, re flmd = l/v = ( dn/dn)/n §. The 
neglect of the coupling between strain and defect den¬ 
sity qualitatively fails in the crystal phases. The full 
compressibility re differs strongly from the approximation 


re°, where both fields are assumed uncorrelated. Thus, 
widely made approximation 39 which identifies re and re° 
fails for cluster crystals. The very close agreement be¬ 
tween re and 1/iq on the other hand, indicates that the 
coupling between strain and density fluctuations is negli¬ 
gible, i.e. n a p « 0; see Sect. IV C 1 Cluster crystals pre¬ 
dominantly accommodate density changes by increasing 
the occupation numbers while keeping the lattice con¬ 
stants fixec£ 17 t With the approximation /i aj g « 0, the 
coefficient ji c a p becomes n° a p ~ v8 a p and the formulas for 
re c and K a simplify to k c w y' 1 + b a $ i ) _1 & r ,5 an d 
re CT ss re ss . Density changes with chemical potential 
similarly at fixed h and er stress tensors. This is in strong 
contrast to the ideal crystal where re 17 equals re c and both 
vanish. 



FIG. 2: The bulk modulus B = - in dimensionless units for 

K 

three different temperatures versus noa 3 . The three points 
are MC simulation result^} 

For a comparison with Monte Carlo (MC) simulations 
the compressibility re from Eq. ( |50| is identified as inverse 
bulk modulus B obtained in Ref~[lB]. Figure[2]shows this 
bulk modulus for three temperatures versus the average 
density. The deviation of the theoretical predictions from 
the simulation data by about 15% is in the same range 
as the deviation of the calculated fcc-bcc transitions from 
the simulated 16 ones; this error is roughly 10%. 

The cluster crystal is an interesting model in the con¬ 
text of defect density fluctuations. The role of the defect 
density is taken by the occupation number n c which ob¬ 
viously is an averaged number; it takes real values, while 
a single lattice site can only be occupied with an integer 
number of particles. There has to be a distribution in 
occupation numbers with the mean value n c and stan¬ 
dard deviation yj < A > which should be connected 
with ( ScSc ). A n c is the occupation number fluctuation 
for each lattice site, so the density <5c(r) has to be inte¬ 
grated over one primitive cell to become equivalent. To 
simplify, we assume that the correlation in occupation 
number/defect density fluctuation vanishes after the first 
Wigner-Seitz-cell, i.e. the occupation number fluctuation 
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FIG. 3: Probability distribution functions for the occupation 
numbers in GEM-4 cluster crystals of fee and bee structure 
from MC simulations^. Ga ussi an distributions with the vari¬ 
ances calculated from Eq. |77|) and the mean value n c ob¬ 


tained through minimization of Eq. (751 (lines) are compared 
with the MC data (symbols). Complete parameters are given 
in table I. 


of each lattice site is independent. With N/n c the num¬ 
ber of lattice sites 


V 


d 3 r(Sc(r)Sc(0 )) 


N 

n c 


(A nl) 


(76) 


This can be rewritten using the compressibility n c from 
Eq.[56| 


(A nl) 



(77) 


Assuming a Gaussian distribution, there is a good match 
of the results for the fee lattice with MC simulationf^, 
as seen in Fig. [3] Table I collects the values of the aver¬ 
ages obtained from the mean-field functional (75) and the 
variances obtained through Eq. Also the percent¬ 

age deviations from the parameters obtained from the 
Gaussian fits to the MC data are shown. The averages 
agree within 1% for both lattices and the variances agree 
for the fee lattice within 10%, which is the same magni¬ 
tude as for the Bulk modulus. For the bcc lattice bigger 
differences between the theoretical and the simulated 17 
occupation number distributions are observed for reasons 
unclear at present. The variances of defect fluctuations in 
bcc and fee crystals are more different in the simulations 
than predicted theoretically. 


C. Dispersion relations and macroscopic density 
correlation function 

The correlation functions for the coarse-grained fields 
can be obtained from the q-dependent elastic coefficients 




V( An c 

) 

n c 


ksT/t 

n oct 3 

MC 

theory 

A[%] 

MC 

theory 

A[%] 

bcc 

1 

6.5 

1.76 

1.32 

33.3 

13.34 

13.24 

0.76 


1.1 

7.5 

1.66 

1.37 

21.2 

15.31 

15.25 

0.39 

fee 

1.1 

8.5 

1.23 

1.31 

6.5 

17.48 

17.49 

0.06 


1 

9 

1.12 

1.24 

10.7 

18.25 

18.44 

1.04 


TABLE I: Variances and averages of the occupation numbers 
in cluster crystals with fee and bcc structure at selected state 
points. The MC results are obtained from Gaussian fits to 
Monte Carlo simulation datcP^ the complete distributions are 
compared in Fig. [3] The theoretical results for the averages 
follow from the mean-field DFT functional (75 I and for the 
variances from 1771). 


according to Eq. |17j). They follow from the density pro¬ 
file obtained by minimizing the DFT free energy func¬ 
tional. The top panel in Fig. [4] shows the dispersion 
relations obtained from diagonalizing the dynamic ma¬ 
trix appearing in the wave-equation of the momentum 
density® D a p(q) = A a g(q) /(mnp) with particle mass 
m, and A given in Eq. ( fl9| ) . A typical state with fee 
lattice is considered. The eigenfrequencies oj of D a p ex¬ 
hibit the familiar longitudinal and (up to two) transversal 
acoustic branches depending on the q-directions followed 
in the first Brillouin zone. Remarkably, the high degree 
of disorder contained in the broad distributions of oc¬ 
cupation numbers does not weaken the solid overly; the 
dispersion relations exhibit the shapes familiar from ideal 
solids and assume magnitudes comparable to the values 
obtained from potential expansions at T = 0 assuming 

idealitjP. 

While the direction-dependence of the dispersion rela¬ 
tions is familiar, the concomitant direction dependence 
of the density correlation functions had not been estab¬ 
lished. The lower panel in Fig. [f] shows the density cor¬ 
relation function (Sn* (q)dn(q)) from Eq. (64) and t%q) 
from Eq. The latter is the q-dependent generaliza¬ 

tion of the thermodynamic derivative v = rip{djji/dn) Uaft 


from Eq. (|39j). Both functions almost completely agree 
for small wavevectors because of the extremely weak cou¬ 
pling between density and strain in cluster crystal s; t he 
coefficient d 2 f /dndu af 3 = fj, a p = [ip5 a p from Eq. (39 d), 
which is diagonal in fee lattices, is very small: [ip/v ss 
2 • 10 -4 . Both functions start deviating for wavevec¬ 
tors approaching the Brillouin zone boundary. Because 
zz(q) possesses a regular small q expansion given in Eq. 
I®, the non-analyticity of the density correlation func¬ 
tion can be brought out by considering the difference 
A(q) = (5n 2 (q))a 3 /V — v~ 1 (q)kBTrip. This A is small 
for small wavevectors because the small factor np enters 
quadratically. Yet, it clearly shows different limits for 
q — > 0 resulting from the direction dependence discussed 
in context with Eq. (62). The insets in Fig. [I] show the 
curves obtained from taking the limit q —> 0 along high- 
symmetry directions in the first Brillouin zone of an fee 
cluster crystal. The directions go from the center T of 
the Brillouin zone along direction [120] (given by Miller 
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FIG. 4: Top panel: Phonon dispersion relation for a cluster 
crystal with fee structure along four symmetry lines in the 
first Brillouin zone; the state at fcsT/e = 1.1 and no<r 3 = 8.5 
is also included in Figs. [2] and [3] and Table I. 

Bottom panel: The q-dependent density correlation func¬ 
tion from Eq. (611 and its dominating part for small q given 
by z/ -1 (q) (blue). Insets: The difference of both quanti¬ 
ties A(q) = {Sn 2 (q))<r 3 /V — u~ 1 (q)kBTnQ in a small range 
around q = 0 for the same symmetry lines. The different 
limits A(q —> 0) depending on direction are apparent. 


indiceiP^) to the point W, along [010] to X, along [111] to 
L, and along [110] to K. Along each of these directions, 
the density correlation function (<5n*(q)(5n(q)) takes a 
different limit for q —> 0. The very small magnitude of 
the differences results from the small value of /io/z' spe¬ 
cific to cluster crystals; the differences are numerically 
reliable. 


D. Discussion of low temperature phase transitions 

Figure [l] shows only little variation of k c iIq, especially 
in the low temperature/high density range. Because 
also varies little, the variance of the occupation number 
fluctuations, (A?r 2 ), is nearly independent of the density 
and scales mainly with the temperature. This points to 
an internal inconsistency of the mean-field description at 
low temperatures. The width of the occupation num¬ 
ber distributions vanishes for T —> 0, yet, non-integer 
average occupation numbers can occur. The failure to 
find integer occupations clearly indicates the break-down 
of mean-field theory for low temperatures. Simulations 
show that the phase diagram of the GEM-4-system ex¬ 
hibits fee phases where the occu patio n numbers take in¬ 
teger values at low temperatures 19 20 . Phase coexistence 
regions lie between them; see Fig. [5] showing simulations 
from Ref. [US]. At critical temperatures each coexistence 
region vanishes, and the homogeneous fee phase with 
a distribution of occupations becomes stable. The MC 
simulation^ suggest that these critical temperatures are 
nearly identical for each phase coexistence, i.e they are 
nearly independent of the density. The mean field den¬ 
sity functional approach only describes the homogeneous 
distributed phase and misses the coexistence regions at 
low temperatures. Potential energy minimization at zero 
temperature gives homogeneous integer occupations and 
rationalizes their coexistence^^. 

Still, the knowledge of the occupation number fluc¬ 
tuations in the homogeneous phase allows to establish 
a criterion when the homogeneous phase is not consis¬ 
tent. We suggest that there is a threshold of the occupa- 



1 1.5 2 2.5 3 

n 0 (J 3 


FIG. 5: Low-temperature phase dia gram of the GEM-4 sys¬ 
tem as determined in MC simulations^* red data points con¬ 
nected by lines as guides to the eye indicate the coexistence 
regions. Pure fee phases with integer site-occupations (de¬ 
noted feen with n = 2, 3,...) survive only at extremely low 
temperatures. Mean-field DFT provides a good estimates of 
the critical temperatures for a reas onable numerical value of 
the occupation number variance, \J (An§) = 0.3 (blue line). 




































17 


tion number variance (A?ij?) where the hopping between 
the lattice sites becomes strong enough to lift the (zero 
temperature) restriction of integer occupation numbers. 
Consequently, for temperatures below this value, we ex¬ 
pect the mean field density functional ( |75[ ) to break down 
and integer occupation phases to become stable, as holds 
at zero temperature. Figure [5] shows that the estimate of 
the occupation number deviation yj (A n%) = 0.3 gives an 
order of magnitude estimate of the critical temperatures. 

The adequacy of the suggested criterion and the sta¬ 
bility of the estimate can be studied in a little more de¬ 
tail. Figure [6] shows the occupation number fluctuation 
for several phases with integer occupations as function 
of temperature. Here integer occupation numbers were 
enforced by hand before minimizing the functional in 
Eq. (751 by varying a only. The critical temperatures 
observed in simulations are well compatible with an oc¬ 
cupation number variation of yj (An;?) ss 0.25, which ap¬ 
pears a rather reasonable value for enabling hopping to 
smear out the occupation numbers on different lattice 
sites. Moreover, varying the value of this criterion moves 
the estimates of the critical temperatures only little. For 
different integer occupations, they differ only slightly. 



ksT/t 


FIG. 6: The occupation number fluctuation yj< A> (stan¬ 
dard deviation) versus the temperature in units of e/fes. The 
standard deviation is a function of temperature and density. 
It is plotted for several integer occupied states. The densities 
are chosen with the approximation n c /no « 2, which differs 
by about two per cent from the optimal DFT-value. The 
red dotted line denotes the point of the curve at the critical 
temperature ksT c /e = 0.471 which is obtained from [N]pT 
simulation^!, the blue line the estimate of Fig. [H] 


VII. CONCLUSIONS AND OUTLOOK 

We derived the thermodynamic expression for the 
isothermal compressibility n in a general crystal, and dis¬ 
cussed its connection to the small wavevector limit of the 
density correlation function. The correlation functions of 


coarse grained fields of macroscopic elasticity theory were 
calculated within the framework of density functional 
theory, allowing for a finite density of defects. Explicit 
expressions for the coefficients in the phenomenological 
free energy in terms of the direct correlation function of 
density functional theory were obtained. The correlation 
function of the coarse-grained density field does not, in 
general, determine the compressibility. For the case of an 
ideal isotropic solid, we could identify the origin of the 
discrepancy from a calculation in macroscopic elasticity 
theory. It arises from the long-ranged strain fluctuations 
which decay like 1/r 3 and thereby cause boundary ef¬ 
fects to enter the elastic energy. While in systems with 
spontaneously broken symmetry, anomalous longitudinal 
correlations exist in general 47 1 (besides the familiar sym¬ 
metry restoring fluctuations^, the present observation 
appears more related to long-ranged dipolar correlations 
in polar fluidiP^. There, the dielectric tensor in response 
to the vacuum electric field depends on the shape of the 
material and on the boundary conditions. It can be con¬ 
nected to a well-defined isotropic dielectric constant only 
via shape/ boundary-effect dependent distribution func¬ 
tions. To work out a corresponding relation for arbitrary 
symmetries and sample shapes of crystalline solids is left 
for future work. 

We applied the theory to the elasticity of cluster crys¬ 
tals made by soft particles. In these crystals, the fluctuat¬ 
ing occupation numbers of lattice sites play the role of lo¬ 
cal defects and strongly affect the stable phases and their 
material responses. Therefore, cluster crystals appear 
an ideal system to test our theory. The obtained com¬ 
pressibilities and occupation number distributions com¬ 
pare well with data obtained in Monte Carlo simulations. 
Mean-field theory breaks down at low temperatures. Yet, 
the theory can be used in order to identify the temper¬ 
ature range where mean-field theory breaks down. This 
provides rather reasonable and stable estimates for the 
critical temperatures, below which the zero-temperature 
phases with integer occupation numbers are stable. 
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Appendix A: Thermodynamic manipulations 


As a consequence of (471 a Gibbs-Duhem relation can 
be derived 


— Vdp + Nd[i + U a pdh a p = 0. (Al) 
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It states that the pressure obeys p = p(p,h a p), which 
can be used to simplify the total differential of the free 
energy density per volume / = F/V = pan — p + u a ph a p. 
It is a proper density because the free energy F is a ho¬ 
mogeneous function of its extensive variables. As a first 
result, from the Gibbs-Duhem relation (Al), the total 
differential of / given in Eq. (491 follows. Also Eq. (Al) 


yields for an isothermal change with dh a p = 0 
Ndp = V dp , 

which can be used for 


(A2) 


1 dN 

A d/7 


V,h aP 


1 dN 
V dp 
_nodV 
V dp 


V.h„ p 


N,h n 


dn 

dp V,h aP 

= non. (A3) 


This verifies Eq. (48) as the compressibility at constant 
stress tensor h a p. 


1. Alternative formula for n and k c 

The discussion of the quadratic terms in the free en¬ 
ergy can be related to more standard thermodynamic 
considerations, which provides additional support for our 
results. To find an alternative formula for the compress¬ 
ibility at constant strain h a p, we start with Eq. (A3) and 
assume a relation u(h, p) 


1 dn 


1 dn 

1 dn 


Tig dp 

h a p 

dp 

u aft ' rig du a p 

v dp 


The last derivative is at constant h a p. With 
dhys 


0 = dh a p = 


du. 


a/3 


du af 3 + dp 


dp 


it can be written as 


du. 


a/3 


dp 


f dhrys 


\du. 


a/3 


-1 


dhrys 


dp 


(A4) 


(A5) 


(A 6 ) 


which leads to the alternative formula 


1 dn 

1 dn 

f 

\ 

1 dh~ys 

rig dp 

u afi rig du af 3 

3 1 \ 


dp 


(A7) 


The inverse of 


is defined bjEl 


d 2 f 


dU a fiUry$ 


= a 


a/ 3^8 


C’ a 0'ysC^g llv — (S a/Jj 6p v + davSpft). 


(AS) 


The unusual definition for the ’’unit matrix” is a con¬ 
sequence from the symmetrisation of the strain tensor, 


u a p = ^{N a up + Vpu a ), and holds for all second order 
derivatives with respect to u a p. 

The thermodynamic derivatives can be expressed 
through the elastic coefficients v, p a p,C a pys as follows: 
The first term of the compressibility is basically the only 
non-vanishing term in the fluid limit 


1 dn 
dp 


/ 2 d P W 1 

v 0 dn U „J 


(A9) 


For the second term the chemical potential p is expressed 
as a function of density and strain tensor p(n,u a p) 


dp 

d d- = 

dn 


dp 
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(A10) 


which yields 
dn 


du. 
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terms 
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— CocpjS l^apn 


dp 


dhys 
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(dp 
ip \dn 


= -non l PyS- 


(A12) 

(A13) 

(A14) 

(A15) 


Now the alternative formula (A7) can be expressed with 


n,p a , 3 , and C%p jg , which yields Eq. (|42). 

The same procedure can be applied to k c 1 which leads 


to 


c 1 dc 

1 dc 

( days 

\ 

1 days 

rig dp 

u afi 1 rig du af 3 

H y dUap 

/V 

dp 


«a/3 

(A16) 

which is an alternative to Eq. (57). The following con¬ 
nection to the elastic constants 


1 dc 
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(A18) 

(A19) 

(A20) 


can be used, to reproduce Eq. (45). 


Appendix B: Elasticity 


With the expressions of Sect. Ill B 2[ the Eq. (63) reads 
(6u*Jup)=Vk B TA-i(q). (Bl) 
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This follows from Eq. (18) with U a (q) = 0, or fi c = 


v8 a p + Tap = 0 in the low q- limit, for a crystal with van¬ 
ishing coupling between strain and defects. If we assume 
the crystal to be ideal, viz. defect free, then addition¬ 
ally the defect density correlations ( ScSc ) and with it k c 
should be zero. This implies that i/(q) _1 vanishes, as 
follows from Eq. (18). The correlations of the coarse¬ 


grained density Eq. (17 1 ) then become 


5n*5n 


VkgTn 2 ^ ~ ( la^ a ^( ( i)Qp — qa(Caef34>qeq<j>) qp, (B2) 


Note, that the limit n c = 0 and Eq. (60) also imply n a = 
0, while k has a finite limit, showing the difference arising 
from the different constant stress tensors. 


The elastostatic theory is contained in the static limit 
of the hydrodynamic equations of motion, see Eqs. (87) 
in Ref. [2]. Without point defects the only non-vanishing 
equation is (87c), which then reads 


where we used the small q expansion of the constants of 
elasticity and took care of the proper symmetric combi¬ 
nation as discussed in [2]. The elastic constants C a g^s o f 
(ideal) elasticity theory correspond to in Eq. (43c I. 

In this ideal crystal approximation, the compressibity be¬ 


comes 


3 

ft - $Ot/3 ocfiryfi) ^ ( (Cij) 


\-l 


i,j =1 


(B3) 


qpC a p-y§qsu^ — 0 . 


(B4) 


But this is just the Fourier-transformed equation of elas- 
tostatics (65). 


1 P.M. Chaikin and T.C. Lubensky, Principles of Condensed 
Matter Physics (Cambridge University Press, Cambridge, 
1995) 

2 C. Walz and M. Fuchs, Phys. Rev. B 81, 134110 (2010) 

3 G. Szamel and M.H. Ernst, Phys. Rev. B 48, 112 (1993) 

4 G. Szamel, J. Stat. Phys. 87, 1067 (1997) 

5 M. Fuchs, in ’Proceedings of the International School of 
Physics ’’Enrico Fermi”’, vol. 184 ’’Physics of Complex 
Colloids”, ed. F. Sciortino, C. Bechinger, and P. Ziherl 
(IOS Press, Amsterdam); also at arXiv: 1209.0389 (2012) 

6 P.C. Martin, O. Parodi, and P.S. Pershan, Phys. Rev. A 
6, 2401 (1972) 

' S. Majaniemi and M. Grant, Phys. Rev. B 75 , 054301 
(2007) 

8 A. Poniewierski and J. Stecki, Mol. Phys. 38, 1931 (1979) 

9 F.H. Stillinger, Jr., Phys.Rev. 142, 237 (1966) 

10 W. Gotze, Phys. Rev. 156, 951 (1967) 

11 W. Gotze and K.H. Michel, Zeit. f. Phys. 217, 170 (1968) 

12 R.F. Kayser, J.B. Hubbard, H.J. Raveche, Phys. Rev. B 
24, 51 (1981) 

13 A. Zippelius, B.I. Halperin, and D.R. Nelson, Phys. Rev. 
B 22, 2514 (1980) 

14 S. Pronk and D.Frenkel, J. Chem. Phys. 120, 67664 (2004) 

15 D. Reinke, H. Stark, H.-H. von Griinberg, A.B. Schofield, 
G. Maret, and U. Gasser, Phys. Rev. Lett. 98, 038301 
(2007) 

18 B. M. Mladek, P. Charbonneau and D. Frenkel, Phys. Rev. 
Lett. 99, 235702 (2007) 

17 B. M. Mladek, D. Gottwald, G. Kahl, M. Neumann, and 
C. N. Likos, Phys. Rev. Lett. 96, 045701 (2006). 

18 Likos et al„ J. Chem Phys. 126, 224502 (2007) 

19 N. B. Wilding, P. Sollich, EPL, 101 10004 (2013) 

20 K. Zhang and P. Charbonneau, Phys. Rev. E 86, 042501 

( 2012 ) 

21 P.D. Fleming and C. Cohen, Phys. Rev. B 13, 500 (1976) 

22 B.J. Berne and R. Pecora, Dynamic Light Scattering , 
(Dover Publications, New York, 2000) 


23 D. Forster, Hydrodynamic Fluctuations, Broken Symme¬ 
try, and Correlation Functions, (Benjamin INC., Reading, 
Massachusetts, 1975) 

24 H. Wagner, Z. Phys. 195, 273 (1966) 

25 J.P. Hansen and I.R. McDonald, Theory of Simple Liquids 
2nd edition, (Academic Press, London, 1996) 

26 J.-L. Barrat and J.-P. Hansen, Basic concepts for simple 
and complex liquids (Cambridge University Press, Cam- 

_ bridge, 2003) 

2 ' J.S. Rowlinson and B. Widom, Molecular Theory of Cap¬ 
illarity, (Clarendon Press, Oxford, 1982) 

28 L.D. Landau and E.M. Lifshitz, Course of Theoretical 
Physics Volume 5, Statistical Physics, (Pergamon Press, 
Oxford, 1970) 

29 J. S. McCarley and N. W. Ashcroft, Phys. Rev. E 55, 4990 
(1997) 

30 F.R. Gantmacher, Applications of the theory of matrices, 
(Interscience Publ., New York, 1959) 

31 R. Evans, Adv. Phys. 28, 143 (1979) 

32 R. Lovett, C.Y. Mou, and F.P. Buff, J. Chem. Phys. 65, 
570 (1976) 

33 M.S. Wertheim, J. Chem. Phys. 65, 2377 (1976) 

34 A.J. Masters, Mol. Phys. 99, 907 (2001) 

35 P. Schofield and J.R. Henderson, Proc. R. Soc. London 

379, 231 (1982); 

J.R. Henderson and P. Schofield, Proc. R. Soc. London 

380, 211 (1982) 

36 N.W. Ashcroft and N.D. Mermin, Solid State Physics, 
(Saunders College, Philadelphia, 1976) 

37 D. Squire, A. Holt, and W. Hoover, Physica42, 388 (1969) 

38 J.F. Lutsko, J. Appl. Phys. 65, 2991 (1989) 

39 D.C. Wallace, in Solid State Physics 25 301, edited by H. 
Ehrenreich, F. Seitz, and D. Turnbull (Academic Press, 
New York, 1970); 

D.C. Wallace, Thermodynamics of Crystals, (Dover, New 
York, 1998) 

40 E. Kroner, Zeit. f. Phys. 136, 402 (1953) 















20 


41 P. De and R.A. Pelcovits, Phys. Rev. B 35, 8609 (1987) 

42 K. Franzrahe, P. Nielaba and S. Sengupta, Phys. Rev. E 
82, 016112 (2010) 

43 J.D. Eshelby, Proc. R. Soc. London A 241, 376 (1957) 

44 K. Garikipati, M. Falk, M. Bouville, B. Puchala, and H. 
Narayanan, J. Mech. Phys. Solids 54, 1929 (2004) 

45 U.M. Titulaer and J.M. Deutch, J. Chem. Phys. 60, 1502 


(1974) 

46 G. Nienhuis and J.M. Deutch, J. Chem. Phys. 55, 4213 
(1971) 

47 W. Zwerger, Phys. Rev. Lett. 92, 027203 (2004) 

48 J. M. Deutch, Ann. Rev. Phys. Chem. 24, 301 (1973) 



